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The problem of search for nearly periodic gravitational wave sources in the data from 
laser interferometric detectors is discussed using a simple model of the signal. Accura- 
cies of estimation of the parameters and computational requirements to do the search 
are assessed. 



1 Introduction 

Pulsars are one of the primary sources of gravitational-waves that can be observed by 
detectors that are currently under construction Bii. The data analysis involved to do the 
search for such sources implies a very heavy computational costB. Here we analyse this 
problem using a simple model of the gravitational-wave signal from a pulsar. A detailed 
summary of the current understanding of the gravitational-wave pulsar phenomenology 
and an analysis of an efficient data analysis technique based on a more accurate model of 
the signal has recently been presented!. 



2 A simple model of the gravitational- wave signal from a pulsar 

The frequency of the gravitational-wave signal from a pulsar will follow its rotational 
frequency and therefore the signal is expected to be almost monochromatic. However the 
amplitude of the signal will be very small and to extract it from the noise we may require 
to integrate the data for several months. Consequently the modulation of the signal due 
to the motion of the detector relative to the solar system barycenter and even very small 
change of the frequency of the pulsar will need to be taken into account. 

Here we consider a simple model of the signal where we take into account only the 
modulation of the signal due to the motion of the Earth around the Sun and we approx- 
imate the change of frequency during the observation time by a Taylor series §. Let Rq 
be 1 astronomical unit (AU), Q = 27r/lyear and let the position of the pulsar on the sky 
be {6, 0) in the coordinate system based on the ecliptic (i.e., 6 = it /2 is Earth-Sun plane, 
and = is position of Earth at t = 0). Let uj{t) be the angular gravitational- wave 
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frequency from the pulsar. We approximate the phase uj{t') dt' of the signal by a power 
series Uit + 0^2^^ + oj^t^ + u^f^ + 0o, where /i = co'i/27r is the frequency of the pulsar at an 
arbitrarily chosen instant of time to and uoi+s is the sth spin-down parameter proportional 
to the sth derivative of the frequency at time to and (po is a constant phase. The number of 
terms needed in the expansion depends on the observation time and the expected values 
of the frequency derivatives. 

We can introduce the following estimates for the spin-down parameters: 

l^l+.l - —X , 1) 

where r is the age of the pulsar and < 1. One can expect that for young pulsars 
are of the order of 1 and less than 1 for old pulsars. Thus we have the following model of 
the gravitational- wave signal: 

h{t) = ho sin[ci;it + uj2t'^ + uJst^ + uj4-^ + uJitQ sin 6 cos(i7t — </)) + (/)(,], (2) 

where ho is the constant amplitude and ^0 = Rq/c. The amplitude ho is estimated asl 

where Izz is the moment of inertia of the pulsar about its rotation axis, r is the distance, 
/i is the gravitational wave frequency and 6 is the ellipticity of the pulsar. The ellipticity 
of 10~^ is an estimate corresponding to the maximum strain that the neutron star crust 
may support. In a realistic model a number of other corrections will need to be taken 
into account!. 



3 Data analysis technique 

The signal given by Eq.(^ will be buried in the noise of the detector. Thus we are 
faced with the problem of detecting the signal and estimating its parameters. A standard 
method is the method of maximum likelihood detection which consists of maximizing the 
likelihood function A with respect to the parameters of the signal. If the maximum of A 
exceeds a certain threshold calculated from the false alarm probability that we can afford 
we say that the signal is detected. The values of the parameters that maximize A are 
said to be maximum likelihood estimators of the parameters of the signal. The magnitude 
of the maximum of A determines the probability of detection of the signal. We assume 
that the noise n in the detector is an additive, stationary, Gaussian, zero-mean random 
process. Then the log likelihood function has the form 

\ogA={x\h)~^-{h\h), (4) 

where x are the data and h is the signal and the scalar product is defined as 

r^™,//, (5) 

•^0 ShU) 

where ~ denotes Fourier transform, * is complex conjugation, and Sh{f) is the spectral 
density of the noise. We can assume that during the time of observation the signal from 
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Table 1: Signal-to- noise ratios for pulsar signals. 





INITIAL 


ADVANCED 


GEO600 


TAMA 


d 


88 


320 


14 
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the pulsar is almost monochromatic. Hence we can approximate the scalar product as 

(^|/^)-^77T / , x{t)h{t)dt, (6) 



where T is the observation time. We can write our signal as h = ho cos 0o^c + sin (pohg- 
Since during the observation time the signal will have many cycles and frequency will not 
change appreciably, to a very good approximation we have 

ih\hs) = 0, (7) 
= {hs\K) = Ho, (8) 

where is a constant ~ Sh{fi) • ^^"^ closed analytic expressions for the maximum 
likelihood estimators of the amplitude ho and the phase 0o of the signal. Substituting 
these expressions into log A and using Eqs.(0) and we get the following formula for 
the reduced likelihood function JF which now depends only on the parameters uji and the 
parameters 9, determining the position of the source in the sky: 

^= ^ ' ' ' • (9) 

For pulsar signal case this last expression can be approximated as 

2 rT/2 fT/2 

Sh{h)T J-T/2 J-T/2 

2 /"^/^ 2 
= ^TTv^l / x{t)exp[-t^^it)-tuit\dt\^ = ———\y\^, (11) 

Oh[jl)J J-T/2 bhijl)^ 

where ^m(t) = uj2t^ + uj^t^ + oj^t^ + uiit^smO cosiVLt — (f)), y{t) is the data multiplied 
by exp[— i$m(^)] and the rectangular window function which is equal to 1 over the time 
interval [—T/2, T/2] and zero otherwise. Tilde denotes the Fourier transform. The above 
calculation suggests one way of evaluating the optimum statistics JF: multiply the data 
by exp[— ^^^(t)] and perform the Fourier transform. This leads to an efficient algorithm 
since we can use the fast Fourier transform. 

The probability of detection of the signal is determined by the signal-to-noise ratio 
d given hy d = (/i|/i)^/^ ~ 's^j?)' "'"^ Table 1 we summarize the numerical values 
for signal-to-noise ratios that can be achieved by laser interferometers currently under 
construction by LIGO, VIRGO, GEO600, and TAMA projects. We choose pulsar with 
Izz = lO^^g cm^, 5 = 10^^. The gravitational-wave frequency /i is 215Hz and the ob- 
servation time T is lO^s. INITIAL assumes approximate model for the noise of LIGO 
and VIRGO detectors at the beginning of their operation. ADVANCED assumes their 
ultimate sensitivity. The distance to the pulsar is taken to be Ikpc except for TAMA 
where it is O.lkpc. 
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The rms errors of the estimators of the parameters of the signal are approximately 
given by the square roots of the diagonal elements of the inverse of the Fisher information 
matrix Fjj given by 

where h' is the signal in terms of parameters 9[ and in our case 

1 /■'^/2 

G--/ cosm,ei)-^t,9',)]dt. (13) 

1 J-T/2 

is called the covariance matrix and it is denoted by C. Instead of the angles 9 and 
it is convenient to introduce the following parameters a and h 

a — u; it Q sin 9 cos (f), (14) 
b — u!itQsm9 sm(f). (15) 

In this new parametrization the phase of the signal has the form 

$(i) = uit + uj2t^ + uj^f + oJif^ + a cos(Ot) + h sin(Oi) + 0^ (16) 

and it is a linear function of the parameters. As a result the correlation function G 
depends only on the difference between the values of the parameters and not on their 
absolute values. Consequently the components of the F matrix are independent of the 
values of the parameters. 

4 Numerical values of the rms errors of the estimators of the peirameters of 
the pulsEir signal 

For observation times T less than about 1/3 of a year one can approximate the components 
of the covariance matrix to a very good accuracy by its leading term in the series expansion 
in T. Let gq^ be the square root of the component Cq^q^ of the covariance matrix. It 
is convenient to express the errors in the parameters uji by the following dimensionless 
quantities 

= ^.7.9xlO-(lf!^)=&(H), (17) 

^1 J- Ji d 

= ^.«.3x:o-(ifi«l^)(j^«12), (1, 
. ^.e,x.o^,ifl)«(l|^,(j^,3(12). (.0, 

The above equations give lower bounds on the rms errors of the spin down parameters. 
The rms error dO, in the position of the source in the sky is given by 

,^ „ f;, sin2f/) , ,l/3yr,i9, IkHz.r, ,10,9 / x 

dQ = 7rsm9aea^ ~ 1.3 x 10"^ ^ Htf^ ^ -7 21 

cos 9 T ji d 

where ag and a^f, are rms errors in the position angles 9 and respectively. 
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5 Computational requirements 



To detect the signal and find the estimators of the parameters we need to find the max- 
imum of the functional T with respect to the parameters. The computational burden 
of the search over the parameter cui is minimized because we can take advantage of the 
speed of the FFT algorithm. The search over the other parameters can be performed by 
means of a bank of filters (templates). The filtered noise in\K) can be thought of as a 
multi-dimensional random process M(^j) with correlation function given by Eq.(|T^). In 
our simple model the correlation function depends only on the difference between the pa- 
rameters and the random process is a generalization of a stationary random process. We 
can generalize the concept of correlation time to such processes, defining the correlation 
hyperellipsoid of the process. The number of independent samples N of such an process 
can be defined as the ratio of the volume of the parameter space V over the area of the 
correlation hyperellipsoid. 

N = r ' (22) 

where n is dimension of the parameter space and Cij is the covariance matrix. 

We use the above formula to estimate the number of independent filters needed to 
probe the signal parameter space. 

For our case since the phase can be eliminated from the search and since to estimate 
the parameter io\ we use the FFT we insert in the above formula the reduced covariance 
matrix which is an n by n submatrix of the covariance matrix corresponding to the n 
parameters that we search for. 

The volume V of the parameter space is given by 

V ^ 7r(^l„,.^e)2u;La.(w)-^(^+')/^ (23) 

where io\max is the maximum frequency we search for and Tmin is the minimum spin-down 
time. 

For observation times less than about 1/3 of a year the number of templates can well 
be approximated by the leading terms of the Taylor expansion of Eq. (^) . 
We obtain the following formulae 

N, ^ 3.4 X 10«(-^)"((!=)i'(i5Z!;). (25) 

l/3yr IkHz r™„ 

iV2 ^ 3.4 X 10l9(-^)l^({^)^(^)^ (26) 

iVs ^ 5.5 X 10^9(-^)2°((^)5(^)^ (27) 
l/3yr IkHz Tmin 

The indicies 1,2,3 mean that 1, 2, and 3 spin-down parameter were included in the 
calculation. The exact values are plotted in Figure 1. 

From Figure 1 we see that at certain observation times the curves intersect. The 
intersection points give times of observation at which we should include a next spin-down 

5 



f_max=lkHz, tau_min=4 0yr 



m 






cu 
+j 




1 9 


rT3 


1 


. 10 


a 




16 




1 


1 U 


0) 




+J 




13 


ll 1 


1 


. 10 


o 




10 


M 


1 


. 10 


0) 




7 


X! 








1 . 10 






s 




10000 . 





































- — — 




_ - - ■ 




^ ^ 




























/ / 












































20 40 60 80 100 120 140 
Observation time [days] 

Figure 1: Number of templates. 



parameter in the search H. To find the number of templates for a given observation time 
one takes the largest of the numbers given above. 

For directed searches for a pulsar of known position in the sky we obtain 

^ ^2. /l 
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The latter formulae are exact. The number of fioating point operations per second (fiops) 
required to perform a search can be obtained by multiplying the above formulae by number 
of operations required to calculate the modulus of the Fourier transform (3A^(log A^+1/2), 
where = 2fimaxT ■, is the number of points of each FFT) and dividing by the time of 
analysis. Assuming that the computation should proceed at the rate of data acquisition 
and for T = 30 days, fimax = IkHz, Tram = 40yr the computing power required is around 
4 X 10^ Tfiops. 
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